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Large-scale Monte Carlo simulations are employed to study phase transitions in the three- 
dimensional compact abelian Higgs model in adjoint representations of the matter field, labelled 
by an integer q, for q — 2,3, 4, 5. We also study various limiting cases of the model, such as the 
Z q lattice gauge theory, dual to the 3DZ q spin model, and the 3DXY spin model which is dual to 
the Z q lattice gauge theory in the limit q — > oo. In addition, for benchmark purposes, we study the 
square lattice 8-vertex model, which is exactly solvable and features non-universal critical exponents. 
We have computed the first, second, and third moments of the action to locate the phase-transition 
of the compact abelian Higgs model in the parameter space (/3, k), where (3 is the coupling constant 
of the matter term, and k is the coupling constant of the gauge term. We have found that for 
q = 3, the three-dimensional compact abelian Higgs model has a phase-transition line j3 c {n) which 
is first order for k below a finite tricritical value K tr i, and second order above. The (5 = oo first 
order phase transition persists for finite f3 and joins the second order phase transition at a tricritical 
point (/3tri,Ktri) = (1.23 ± 0.03, 1.73 ± 0.03). For all other integer q > 2 we have considered, the 
entire phase-transition line (3 c (k) is critical. We have used finite-size scaling of the second and third 
moments of the action to extract critical exponents a and v without invoking hyperscaling, for the 
XY model, the Z2 spin and lattice gauge models, as well as the compact abelian Higgs model for 
q = 2 and q = 3. In all cases, we have found that for practical system sizes, the third moment gives 
scaling of superior quality compared to the second moment. We have also computed the exponent 
ratio for the q = 2 compact U(l) Higgs model along the critical line, finding a continuously varying 
ratio (1 + ct)/v, as well as continuously varying a and v as k is increased from 0.76 to 00, with the 
Ising universality class (1 + a)/u — 1.763 as a limiting case for (3 — > 00, k — > 0.761, and the XY 
universality class (1 + a)jv = 1.467 as a limiting case for (5 — > 0.454, ft — * 00. However, the critical 
line exhibits a remarkable resilience of Z2 criticality as f3 is reduced along the critical line. Thus, the 
three-dimensional compact abelian Higgs model for q — 2 appears to represent a fixed-line theory 
defining a new universality class. We relate these results to a recent microscopic description of 
zero-temperature quantum phase transitions within insulating phases of strongly correlated systems 
in two spatial dimensions, proposing the above to be the universality class of the zero-temperature 
quantum phase transition from a Mott-Hubbard insulator to a charge-fractionalized insulator in two 
spatial dimensions, which thus is that of the 3-D Ising model for a considerable range of parameters. 

PACS numbers: 74.20.-z, 05.10Cc, 11.25. Hf 



I. INTRODUCTION 

Lattice gauge theories in 2+1 dimensions with compact 
gauge-fields coupled to matter fields, have recently come 
under close scrutiny as effective theories of strongly cor- 
related fermion systems in two spatial dimensions at zero 
temperature. Phase transitions in such three dimensional 
models then correspond to quantum phase transitions in 
a system at zero temperature in two spatial dimensions. 
A central issue is whether such systems of strongly corre- 
lated fermions can suffer quantum phase transitions from 
Fermi-liquid metallic states to states where the quasipar- 
ticle concept has broken down and given way to singular 
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Fermi liquids or electron-splintered states 0, Q . Such 
quantum phase transitions may be related to phase tran- 
sitions such as confinement-deconfinement transitions in 
2 + 1 dimensional compact gauge theories. This fact has 
resulted in focused attention on effective gauge theories 
of matter fields representing charge doped into Mott- 
Hubbard insulators, coupled to fluctuating gauge fields 
representing strong constraints on the dynamics of the 
fermions on the underlying lattice on which the models 
are defined 0, HI IS 0> HI • Since these are lattice models, 
the corresponding gauge fields are necessarily compact. 
Compact U(l) gauge fields in 2 + 1 dimensions support 
stable topological defects in the form of monopole con- 
figurations, and it has been suggested that the unbinding 
of such monopoles, a confinement-deconfinement transi- 
tion, may be relevant for such phenomena as spin-charge 
separation in strongly correlated systems [S IS El ■ Con- 
finement here refers to the confinement of test charges 
in the problem, not of topological defects of the gauge 
field (which are space-time instantons, and will hereafter 
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be referred to as "monopole" configurations). Alterna- 
tive formulations in terms of a lattice Ising gauge theory 
[l?M c oupled to matter fields, have also been put forth 
[21 llll . This is largely motivated, it would appear, by 
the fact that it is highly controversial whether a 2 + 1- 
dimensional U(l) gauge theory with matter fields in the 
fundamental representation will undergo a confinement - 
deconfinement transition. In the absence of matter fields, 
compact U(l) gauge theories are known to be perma- 
nently confined in d = 2 + 1, while the pure Z 2 gauge 
theory is known to have a second order phase transition in 
the inverted 31? Ising universality class. It is sometimes 
stated, without much justification, that the presence of 
matter fields in the fundamental representation will not 
alter the picture that emerges in their absence, both for 
the U(l) case and the Z2 case in 2 + 1 dimensions. 

One major problem that arises in this context, is that 
the Wilson loop non-local gauge-invariant order param- 
eter, which has proven itself to be very useful in the ab- 
sence of dynamical matter-fields to distinguish confined 
from deconfined phases, is rendered useless by the pres- 
ence of them. In particular, when the symmetry group of 
the matter field is contained in the symmetry group of the 
gauge field, one can demonstrate rigorously under oth- 
erwise quite general conditions, that the Wilson loop is 
bounded from below by a perimeter law, not an area law, 
despite the fact that models with this property definitely 
has phase transitions from confined to deconfined phases 
[l3| . Hence, the Wilson loop and the related Polyakov 
loop, are no longer useful order parameters for the prob- 
lem. Physically, this is due to the fact that in the pres- 
ence of a dynamically fluctuating matter-field coupled to 
the gauge field, particle-hole excitations are generated 
from the vacuum and will effectively screen the interac- 
tion between two static test-charges introduced into the 
system, which the Wilson loop is a measure of. Hence, a 
perimeter law is always obtained. 

Recently, two of us have shown that a confinemcnt- 
deconfinement transition may take place in three- 
dimensional compact U(l) gauge theory coupled to mat- 
ter fie lds, when the matter field exhibits critical fluctua- 
tions [Tj, [l5j . Such matter-field fluctuations endow the 
gauge-field with an anomalous scaling dimension r/A with 
a value (l6lll7| . 

VA = 4-d, (1) 

where d is the dimensionality of the system. This value 
for t]a is protected by gauge-invariance. In d = 3, this 
transforms the gauge-field propagator in a striking man- 
ner such as to allow a confinement-deconfinement tran- 
sition to take place via a three-dimensional Kosterlitz- 
Thouless like phase-transition of unbinding of point-like 
monopole configurations of the gauge-field |Tj, UJ| ■ The 
treatment of Refs. closely parallels that of Ref. 
ITsI and is based on a dual description in the continuum. 
However, it is far from obvious what the corresponding 
lattice gauge theory, if any, that would yield such results, 
could be. 



In the absence of a clear-cut order-parameter criterion 
for distinguishing various phases of such matter-coupled 
gauge fields, it would be advantageous to be able to dis- 
tinguish various phases by direct "thermodynamic" mea- 
surements reliably exhibiting possible non-analytic be- 
havior. In this paper, we will use one such measurement 
recently introduced by us |l9j . namely finite-size scaling 
of the third moment of the action of the lattice model. 
This will turn out to be a superior quantity to study for 
this purpose, compared to the second moment. It brings 
out non-analytic thermodynamics and precision values of 
the specific heat critical exponent a and the correlation 
length exponent v through finite-size scaling analysis per- 
formed on practical system sizes . The third moment 
has the advantage of not being contaminated by contri- 
butions from the second moment. The second moment 
is known to be a notoriously difficult quantity to use in 
particular for extracting specific heat exponents, due to 
large corrections to scaling coming from confluent sin- 
gularities for practical system sizes. Moreover, the third 
moment is the simplest quantity to compute which has an 
extra feature which even moments do not have. It has a 
double-peak structure where the width between the peaks 
also exhibits scaling. This allows us to extract separaely 
two exponents, a and u, from measurements of the third 
moment alone, without having to invoke hyperscaling. 

A lattice model of particular interest in this context is 
the abelian U(l) Higgs model with a compact gauge-field 
[3 coupled minimally to a U(l) bosonic matter field [23, 
|21| with a gauge-charge q. It is defined by the partition 
function given by the following functional integral 




exp [ft H$ + k H K ] 



H p = ^(1 - cos(0 3> )) 

H K = 5^(1 -cos^)), (2) 

where N is the number of lattice sites, and we have de- 
fined 

Sj> = A^j - qA jti 

— £^i>\AvAj\. (3) 

Here, e^ u \ is the completely antisymmetric tensor. More- 
over, V. , denotes a sum over sites of the lattice, while 
— jit 1 

J2p n denotes a sum over the plaquettes of the lattice. 
We will use the variables (x = 1/[k + 1], y = 1/\B + 1]) 
when discussing the possible phases of this model |2(j . In 
Eq. (J2J, 6 is the phase of a scalar matter field with unit 
norm representing holons, is a forward lattice differ- 
ence operator in direction fi, while A, M is the fluctuating 
gauge field enforcing the onsite constraints reflecting the 
strong correlations in the problem. We are neglecting 
amplitude fluctuations of the matter fields, working in 
the "London" -limit. 
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Let us summarize what is known about this model. 
When q = 0, the matter field decouples from the gauge 
field. It is well known that the model then has one critical 
point in the universality class of the 3DXY model, and 
in the Villain approximati on ( which is most often used 
when dualizing the model [2(l|). the critical value y c is 
given by y c w 3/4 [22l l23l |. On the other hand, the pure 
gau ge theory is permanently confined for all values of k 
[l8| . Consider next q = 1. Then, Eq. (|5J) is trivial on the 
line x = 1,0 < y < 1 with no phase transition for any 
value of y. On the line < x < 1, y = 1 the matter field 
is absent and the theory is permanently confined |18| . 
For a further enumeration on rigourous results both on 
the non-compact and compact version of this model, see 
also 0. 

For arbitrary {[3, k), the case d = 3, q = 1 is con- 
troversial. It is however clear that no ordinary second 
order phase transition with a local order parameter ex- 
ists for the model in this case. When matter-fields are 
coupled to a compact gauge-field in a continuum the- 
ory, the permanent confinement of the pure gauge theory 
is destroyed and a confinement-deconfinemcnt transition 
may take place via a Kosterlitz-Thouless like unbinding 
of monopole configurations [l4| in three dimensions. This 
is due to the appearance of an anomalous scaling dimen- 
sion of the gau ge-fi eld induced by critical matter-field 
fluctuations [la. \u\ . The role of the anomalous dimen- 
sion has also been studied recently at finite temperature, 
in pure compact QED in d = 3 with no matter fields 
present |2^. In both Refs. and l25t the anomalous 
scaling dimension plays a crucial hole. However, in Ref. 
I2H the deconfinement transition occurs due to finite tem- 
perature, and there are no matter fields. 

The role of an anomalous scaling dimensions in driv- 
ing a recombination of monopole defects of the gauge 
field into dipole configurations, corresponding to a 
confinement-deconfinement transition, has been studied 
numerically for the case q = 1 in the presence of matter- 
fields [26J . The authors of Ref. |26| in this case reach con- 
clusions in agreement with Ref. Thus, even putting 
the issue of whether this signals a phase-transition or not 
aside, it is clear that also matter fields in the fundamen- 
tal representation do matter in this problem. It is also 
interesting to note in this context, that some time ago a 
rather remarkable paper |27| , a Kosterlitz-Thouless tran- 
sition was claimed in a three-dimensional theory of inte- 
ger point charges interacting via a logarithmic potential. 
There however, the origin of the logarithmic interaction 
was due to higher order anisotropic gradient terms, es- 
sentially an input to the theory, and not a result of an 
anomalous scaling dimension appearing due to criticality 
in the mass-sector of a compact U(l) gauge theory. This 
however, contrasts with the case considered in Ref. IhH 
The anisotropy effectively leads to a dimensional reduc- 
tion and a resulting standard Kosterlitz-Thouless phase 
transition in two dimensions. 

For q > 1 the model in Eq. J2Jl exhibits a behavior that 
completely sets it apart from the case q = 1. The reason 



is that for q > 1, the matter field is no longer in the 
fundamental representation. Thus, in contrast to the no- 
compact case of the standard Maxwell term of the gauge 
sector, the gauge charge cannot be simply scaled away. 
The case q — 2 is particularly interesting in the context of 
electron fractionalization [23] ■ This is a case that we will 
consider in detail in this paper. This theory arises as a 
special limit of a bosonic model exhibiting fractionalized 
phases considered recently [2^. The phase diagram for 
q = 2 was briefly discussed long ago by Fradkin and 
Shenker [2l|. For d = 3 the phase diagram is divided 
into two phases, a confined and a deconfined-Higgs phase. 
There is no Coulomb phase in three dimensions. The 
model becomes a Z2 gauge theory on the line y = 0. 
This limit suffices to bring out the fundamental difference 
between the q = 1 and q — 2 cases, since on the line 
y = 0, the q = 1 case is trivial. We will discuss this 
in more detail below. The vortex content of the q = 2 
case is different: Z2 vortices, or visons, may arise in the 
deconfined-Higgs phase. Due to the visons the flux is 
quantized in units of 27r instead of 7r. This means that 
the excitations in the deconfined-Higgs phase have charge 
e/2. In the context of the bosonic model in Ref. l28l 
this phase corresponds to a fractionalized insulator. The 
confined phase, on the other hand, features excitations 
with charge e and should correspond to a conventional 
Mott insulator [28[. Thus, the model J5J) can be thought 
as describing a particular case of insulator -fractionalized 
insulator transition. The full bosonic model considered in 
Ref. |28| has a more rich phase structure depending on the 
values of the parameters. For example, a superfluid phase 
is also possible. We will not consider such a situation 
here, but only a special limit of the bosonic model of 
Ref|2l 

The purpose of this paper is to make a detailed nu- 
merical study of the phase structure of An impor- 
tant point concerns the universality class of the phase 
transition. This will be the main topic discussed in this 
paper. We have performed a large scale Monte Carlo 
study which gives a very complete picture of the univer- 
sality class of the transition. The present study com- 
plements and goes far beyond our previous lar ge s cale 
Monte Carlo study of the model |19|. In Ref. we 
have shown that the Ising [Z2) universality class domi- 
nates over a significant portion of the critical line of the 
theory. Moreover, we have shown that there is a region 
of parameters where the critical exponents are continu- 
ously varying. This was interpreted as evidence for the 
existence of a line of fixed points in the renormalization 
group (RG) flow diagram. To the portions of the critical 
line corresponding respectively to Z2 and XY critical be- 
havior, are associated Ising and XY fixed points. Besides 
these two fixed points there exists a fixed line correspond- 
ing to a critical phase. As discussed recently 0, electron 
fractionalization in 1 + 1 dimensions is associated to a 
quantum critical phase. Ordinarily, in more than one 
spatial dimension quantum phase transitions are associ- 
ated with a critical point. Thus, we would expect that 
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electron fractionalization would occur at a critical point. 
Our analysis clearly shows that a critical phase exists for 
the model at q = 2. This becomes particularly clear 
on the portion of the critical line corresponding to large 
values of n, where we find that the critical exponents 
vary smoothly with coupling constants, approaching the 
3DXY value only slowly. The three-dimensional KT-likc 
scenario [lj, [l5j for the q = 1 model is also an example 
of critical phase occuring in higher dimensions. However, 
due to the vortex content of the model at q = 1, it is not 
entirely obvious that the corresponding deconfinement 
transition is really associated with electron fractionaliza- 
tion 0. 

In this paper the cases q — 3,4, and 5 will also be 
considered. The case q — 3 is particularly interesting. 
While for q = 2, 4, and 5 the entire line separating the 
two phases are critical, this is not true when q = 3. In 
this case there is a value of k below which the transition 
is first-order, being second-order otherwise. This point 
where the transition changes from second order to first 
order, is clearly a tricritical point. We emphasize that 
among the situations analysed in this paper, only the case 
q = 3 exhibits a tricritical point. When q = 3, charge is 
fractionalized in such a way that excitations carry charge 
e/3. This situation is reminiscent of the v = 1/3 state 
in the fractional quantum Hall effect jsij. In general 
the charge for arbitrary q will fractionalize as e/q. This 
gives an elementary flux quantum cf>o — 27rq/e, contrast- 
ing with the more familiar situation of flux quantization 
in a superconductor. Indeed, if we take the example of 
the q = 2 theory, we see that the flux quantum is dou- 
bled, while in a superconductor it is halved because the 
Cooper pair has charge 2e. Other interesting possibil- 
ities of vortex/charge fractionalization were considered 
recently in certain Josephson junction arrays [3lj . 

In all these examples, it is readily seen that the type of 
flux quantization involved depends on specifying whether 
charge is fractionalized or not. Therefore, in order to de- 
tect charge fractionalization, experimental measurements 
of flux as the one proposed in Ref. 32 in the context 
of underdoped cuprates are important. The experiment 
proposed in Ref. l33has now been performed [3j|, but has 
failed in detecting electron fractionalization. The failure 
of this experiment does not, however, entirely rule out 
charge fractionalization in the cuprates for the following 
deep reason. In the zero doping insulating regime, the 
system can effectively be described by the antiferromag- 
netic quantum Heisenberg model. This model is known 
to have a local SU(2) gauge symmetry [24]]. Upon dop- 
ing, this symmetry is broken down to U(l). Since this 
U(l) is a subgroup of SU(2), it is compact and carries 
the "charge" label q corresponding to the different rep- 
resentations. Within the q — 2 represention, the double 
broken symmetry pattern SU(2) — > U{\) — > Z<i is possi- 
ble. In such a situation, there would be no stable vison 
and any experimental attempt based on a vison trapping 



experiment would fail. However, this does not mean that 
charge is not fractionalized. The Z2 residual subgroup 
will imply the existence of magnetic monopoles which, by 
duality, accounts for the existence of fractional charges 
e/2. 

In the perspective of such a possible scenario, alterna- 
tive experimental proposals should be envisioned and the 
identification of the universality class of possible phase 
transitions the q = 2 case is important in this context. 

The plan of the paper is the following. In Section II, 
we discuss the limits of the model, pointing out the dif- 
ferences with respect to the q — 1 case. In Section III, 
we explain the finite-size scaling analysis of the third mo- 
ment of the free energy used in this paper. In Section IV 
we discuss the details of our large-scale Monte Carlo sim- 
ulations. As a check on the quality of the simulations, 
we have performed several benchmark simulations which 
we present in parallel with our new results on the model 
defined by Eq. J5J). Section V concludes the paper. 
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LIMITS OF THE COMPACT GAUGE 
THEORY 
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00, k finite 



In terms of the formulation of the model in Eq. J3J, 
this limit leads to the constraint A^Oj — qAj^ = 2irli^ 
where li^ is integer valued. Substituting this into the 
gauge-field term, we find 



A' 



n e ex p 

j = l (j,n = -oo 



>«^(l-cos(— £ J> )) 

P,M q 



(4) 



where £ JiAt = e^xA^lj^x is also integer valued. For q = 1 
the model is therefore trivial in this particular limit, in 
accordance with the discussion in the previous section. 
For q = 2, the model is equivalent to the lattice Z2 gauge 
theory [l(J, and the critical point of the model in this 
limit is thus in the inverted Ising universality class (in 
analogy with the inverted XY universality class of the 
dualized 3DXY model [13 )■ Therefore a and v are those 
of the global 3D Ising model. For integer q > 3, the 
critical exponents a and v we will find in the limit (3 — > 
00, k fixed, will be those of the Z q spin model, since as 
we shall see below, the Z q gauge theory is dual to the Z q 
spin model in d = 3. 

For arbitrary integer-valued gauge-charge q (i.e. la- 
belling of matter-field representation where q = 1 means 
fundamental representation, while q > 2 means higher 
representations), we may write the action in Eq. J5J in 
the Villain approximation, replacing the cosine-terms by 
periodic quadratic parts, after which the model may be 
written in terms of its topological defects as follows |a.l20| 
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Z = Z E E S ^uJju,9Qj ex P 



3,k ^ 



" t -^QjQk }D(j ~ k,m 2 ) 



(5) 



where D(j—k, m 2 ) — (— A 2 x +m 2 )~ 1 6jk andm 2 = q 2 (3/n. 
Here Zq is the partition function for massive spin waves 
ppf. and is an analytic function of coupling constants 
which will be omitted from now on. Note the appearance 
of the constraint 

A v J jv = qQ j ; qeN (6) 

in the summation, which will be important in what fol- 
lows. Here, Qj € Z is the monopole charge on (dual) 
lattice site number j, while Jj V are topological currents 
representing segments of open-ended strings terminating 
on monopoles, or closed loops. These are the only stable 
topological objects of the theory [20|. For a recent treat- 
ment of the interplay between abelian monopole conden- 



sation and vortex condensation in lattice gauge theories, 
see [36) . 

In the limit f3 — ► oo at fixed K, the partition function 
in Eq. JSJ) takes the form 

^=E E^^expf-^^^y (7) 
W : :./ • V / 

This is easily seen to be the loop-gas representation of the 
global Z q theory in the Villain approximation [37l | . Using 
the integral representation of the Kronecker delta and 
summing over Jj v using the Poisson formula, we obtain 
up to an overall non-singular factor 



*= E 



2n 



exp 



8tt 2 k 



]T(A^ 



2nN JU ) 2 



(8) 
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By employing the Poisson identity 



J2 J qQe = ^ E w - 27rZ /«)' ( 9 ) 

Q— — QO l— — oo 

we obtain 



Z 



E 



{ir- 



E ex p 



8tt 2 



E 



2tt 



-AJi - 



(10) 
model 



The above is precisely the Villain form of a 

M 

Since Eqs. @ and J2J are dual (up to a Villain ap- 
proximation), and Eq. Q is a loop-gas representation of 
the global Z q theory while Eq. (0} is the Z q lattice gauge 
theory, we conclude that the global and local Z q theories 
are dual to each other in d = 3, in agreement with the 
results of Ref. [3^. Hence, again we conclude that the 
model Eq. J5J) in the limit (3 — > oo, n fixed, should have 
critical exponents a and v consistent with the g-state 
clock model universality class. Since lim^oo Z q = U(l), 
the above fits nicely in with what is known for the U(l) 
case in d = 3, where the global 3DXY model dualizes 
into a U(l) gauge theory |H S E3, EH • 



From Eq. J7J), it is seen that the cases q = 1 and 
q =/= 1 are fundamentally different. For q = 1, the 
summations over {Qi} may be performed to produce a 
unit factor at each of the N dual lattice site, thus com- 
pletely eliminating the constraint j23- Hence, we have 

/ 2 \ N 

Z = I $3(0, e -2 *" K ) ) where the elliptic Jacobi functions 

are given by #3(2, q) = X)^L-oo 1 n exp(27rinz), which 
are analytic functions and hence no phase-transition oc- 
curs at any value of k for q = 1 in this limit. For (integer) 
q > 2, a phase transition is known to survive |2l| . The 
restoration of a phase-transition for integer q > 1 is, in 
this language, crucially dependent on the presence of the 
constraint Eq. Jfjjl for q 7^ 1. Even when one sums over 
all possible values of {Qj}, this still represents a real con- 
straint on the vortex-configurations of the system, since 
it cannot be eliminated by summation as in the case of 
q = 1. This suffices to convert the theory to a strongly 
interacting one capable of sustaining a phase-transition, 
in contrast to the effectively (discrete Gaussian) nonin- 
teracting case q = 1. 
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B. f3 finite, k — > oo 

We now discuss the limit (5 finite, n — > oo. Note from 
Eq. ©, that when n — > oo, we have e^xA^Ajx — 
27tMj jM , where Mj^ is integer valued. This shows that in 
this limit, the gauge-field fluctuations are precisely like 
transverse phase- fluctuations in the 3DAY"-model [42^ . 
and the integer Mj tfl plays the role of vorticity. Thus, 
even when k — > oo, gauge-field fluctuations are not com- 
pletely suppressed, we are not allowed to set Mj ^ = 0, 
and hence we cannot set Aj tfi — A^Xj m the nrs t term in 
Eq. (J2J). However, we may write the partition function 
on the form 



Z 





/" 


TT ^A/> 
J-i 2tt 




J — 7T 





exp[/3ff e ], (11) 



where the functional integral over the gauge-field is now 
constrained. Introducing the Villain-approximation, us- 
ing the Poisson summation formula, integrating over 9j , 
and performing a partial integration, leads to the parti- 
tion function which is equivalent to an ordinary loop gas 
with completely suppressed gauge-field fluctuations when 
q is an integer. The gauge-field fluctuations of the com- 
pact sector produce closed vortex loops with long-range 
interactions such that the loop-tangle is incompressible, 
entirely equivalent to ordinary matter- field vortex loops. 
No screening of their test charges is produced by these 
constrained gauge-field fluctuations and monopole config- 
urations are suppressed. Hence, the model in this limit 
is, after all, equivalent to the 3DXY model, as it also 
is in the non-compact case where gauge-field fluctuations 
are entirely suppressed. 

Ideally, to explore the entire phase diagram of the 
model Eq. (0), we could carry out Monte Carlo (MC) 
simulations in conjunction with finite-size scaling analy- 
sis of standard quantities such as the susceptibility or the 
specific heat. The problem with the former quantity is 
that we have no candidate gauge-invariant order param- 
eter for this model, in view of the results of [l3l l43l l44| . 
The problem with the latter is that the second moment 
of the action is marred with huge corrections to scaling 
and is hence unsuitable to extract critical exponents from 
simulations on practical system sizes. To circumvent this 
difficulty, we will consider various higher order moments 
of the action appearing in the partition function in or- 
der to obtain quantities with good scaling for practical 
system sizes we can handle in Monte Carlo simulations. 



III. MOMENTS OF THE ACTION 

Given a singular contribution to the free energy F = 
In Z of a system with an inverse temperature coupling 



0, 



where h plays the role of a scaling variable, A is 
some scaling exponent, lim^oo ^±(2;) ~ x( 2 ~ a ^ A , 
\im x ^Q J 7 ±(x) = A±, A± are critical amplitudes, and 
t = {/3 — Pc)/(3 is some deviation from a critical coupling, 
then the non-analytic contribution to the susceptibility 
of the action is given by 



c~\t\- a C ± (h/\t\ A ), 



(13) 



where \im x _, 00 C±(x) ~ x a / A , and lim x ^oC±(x) = A±. 

At a critical point, this quantity will scale as C ~ L a ' v , 
where the volume of the system is given by L x L x L, 
provided the system exhibits one diverging length scale 
at the phase-transition. Here, v is the correlation length 
critical exponent, £ ~ \t\~ 1 ' close to a critical point. A 
problem arises if a < 0, as in the 3DXY model, since 
one then gets an increasing peak in C itself with increas- 
ing L, which however eventually will no longer scale with 
L. Unfortunately, but quite typically, impractically large 
system sizes are needed to eventually distinguish correc- 
tions to scaling from actual scaling in the second mo- 
ment, particularly so when a < 0. Thus, C exhibits a 
finite cusp non-analyticity which does scale, superposed 
on a large regular background which eventually will not. 
It would be advantageous to be able to subtract out this 
background, or at the very least bring out the leading 
non-analyticity dominating the scaling more clearly rel- 
ative to confluent singularities, or corrections to scaling. 
This can effectively be achieved by taking one further 
derivative of the action with respect to the coupling con- 
stant as follows 



d 3 F sing 
d(3 3 



\t\-^ g±(h/\t\ A ), 



(14) 



where the scaling function Q±(x) has the following prop- 



erties lim-E^oo G±(x) 



-(l+a)/A 



, and lim-r^o G±(x) 



A±, which will scale as L^ 1+a ^ 1 ' at a critical point. 
More generally, the n-th moment 



r) n F ■ 

u r SI 



d(3 n 



| t |-(n-2+«) Q ± (h/\ t \% 



(15) 



F sins ~\t\ 2 - a F±(h/\t\ A ), 



(12) 



will scale as i(™ _2 + Q )/ !y at a critical point, and therefore 
by computing two moments, say the third (n = 3) and 
the fourth [n = 4), it is also possible to extract a and 
v separately without utilizing the hyperscaling relation 
a — 2 — dv. 

In fact, this may be obtained from the third moment 
M3 alone, since the width between the negative and pos- 
itive peaks scales as L~ x l v . Thus M3 yields independent 
measurements of both (1 + ex) fv and 1/v. The above 
procedures may serve as checks on the validity of hyper- 
scaling. This is known to be violated above the upper 
critical dimension in spin-models and systems with long- 
range interactions due to the presence of dangerous irrel- 
evant variables |45| . Long-range interactions are features 
of the dual models of some of the theories we consider and 
little is known about the presence of dangerous irrelevant 
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operators in some of the models we study in this paper. 
Hence, caution is necessary in extracting individual ex- 
ponents. At any rate, the above seems to be a useful 
way of extracting two exponents a, v from measurements 
of one quantity, namely M3. The procedure described 
above will presumably turn out to be useful in a host of 
other models in statistical physics. 



M, 



K H 




^^^^^ 


£(l+a)/v 



Coupling constant 

FIG. 1: Generic third moment of action, M3, showing how 
finite-size scaling is used to extract a and v. 

Consider next a model with two coupling constants X\ 
and x 2 defined by its action 



S = x\ Hi 
In this case, the scaling of 



x 2 H 2 . (16) 
could conceivably de- 



dxidxj 

pend on the direction in parameter space in which the 
critical line is crossed. To take this complication into ac- 
count, we consider contributions to the free energy to sec- 
ond order in devitations in the coupling constants from 
their critical values, by expanding to second order around 
the critical point 



F = -InZ 

1 S 2 F 

= F({x c })+ 



{*c} 



This contribution to the free energy can be written 
F( 2 \ where Fa is the "fluctuation matrix" 



5x T FijSx = 
defined by 



F, 



(Hi HA - (Hi) 



d 2 F 



(17) 



This can be diagonalized via an orthogonal transfor- 
mation by rotating the original coupling constant basis 
through an angle 9 46j, providing a new set of coordi- 
nates which are linear combinations of coupling constants 



x 1 — x± cost; + X2 sin ft 
x 2 = —xi sin 9 + X2 cos 9, 

with the corresponding uncorrelated operators 

H[ = Ht cos 9 + H 2 sin 9 



HL 



-Hi sin6> + H 2 cos 9. 



(18) 



(19) 



This yields 



(Sx'i,5x' 2 [ 







A_ 



Sx'i 
5x' 



(6x'i) 2 \ + + (6x' 2 ) 2 \-, 



(20) 



where A+ and A_ are the larger and smaller eigenvalue, 
respectively. If A + ^> A_ the first term will dominate 
the leading order corrections to F({x c }). Hence, to ob- 
tain proper scaling, the second derivative of the free en- 
ergy should be evaluated along the direction of the cor- 
responding eigenvector, 



d 2 F 

^75~K-x lc | Q , 

when crossing the critical line. 

From Eq. J5J), we get the standard expressions 

dF 
dpi 
d 2 F 



(21) 



= -<H U 



dpi 2 



(22) 



where p is the rotated coordinate along the eigenvector 
corresponding to the largest eigenvalue, and where we 



have defined G M = 
Hu defined by Eq. 



- (H^ with rotated operator 
The above expressions essen- 



tially represent generalizations of the expressions for in- 
ternal energy and specific heat. The second derivative 
contributes to C M which we will refer to as the singular 
part of the second moment of the action, or the action 
susceptibility. Similarly, we obtain 



d 3 F 



3 — (G>)> 



and in general we have for the n-th moment 
d n F 



(23) 



(24) 



IV. MONTE CARLO SIMULATIONS 
A. Details of the MC simulations 

The critical properties of the models are investigated 
using large scale Monte-Carlo simulations in conjunction 
with finite-size scaling analysis. A Monte-Carlo move is 
an attempt to replace the field value at one particular 
point by a randomly chosen value. The move is rejected 
or accepted according to the standard Metropolis algo- 
rithm. In the case of the XY-model, the phase 9j is the 
relevant field, whereas in the full abelian Higgs model 
both the phase and the three components of the gauge 
field Aj iM are subjected to the Metropolis algorithm. One 
sweep consists of traversing the system L x L x L while 
attempting a Metropolis update of each field component 
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once. The acceptance rate of the Metropolis algorithm is 
kept fixed between 60% and 70% by dynamically adjust- 
ing the maximum allowed changes in the fields. However, 
for the model defined in Eq. J2J at large f3, where Z q 
gauge behaviour is expected, the algorithm for control- 
ling the acceptance rate has been relaxed by fixing the 
maximum allowed change to 2tt. There is no gauge fix- 
ing involved in these simulations, and periodic boundary 
conditions are used in all directions. The excess gauge- 
volume due to the summation in the partition function 
over (redundant) gauge-equivalent field configurations, 
will cancel out when computing all averages of gauge- 
invariant quantities, in particular when computing mo- 
ments of the gauge-invariant action. 

The third moment of the action will typically behave 
as shown in Fig. ^ The MC simulations for the com- 
pact abelian Higgs model will be performed for a set of 
coupling constants that span a line across the phase tran- 
sition. We diagonalize the "fluctuation" matrix Eq. i|17|) 
and simulate along the trajectory (3{n) — (3 C + a(n — k c ), 
with a being determined from the eigenvector with the 
largest eigenvalue. In all models that we consider in this 
paper, including benchmarks models, we then compute 
the moments of the action according to Eq. I|24|l applying 
Ferrenberg-Swendsen multihistogram reweighting analy- 
sis with jackknifing error estimate. For the benchmark 
models, this procedure produces curves as shown in Figs. 
|3]and0] The top and bottom values as well as their posi- 
tions for different system sizes are then used to produce 
scaling plots as shown in Fig. [5j The combinations of 
critical exponents (1 + a)/v and \ jv are then extracted 
by bootstrap regression analysis. 



B. Simulations of the 8-vertex model 

As a first benchmark on our method of extracting the 
exponents v and a separately, we consider the 8-vertex 
model j47j on a square lattice. This model has the virtue 
of being exactly solvable, and hence an analytic expres- 
sion for the exponent v is known. Moreover, it has crit- 
ical exponents that are non-universal [47|. The 8-vertex 
model on a square lattice ma y be formulated as a gener- 
alized Ising model as follows |47| 

Zsv = ^2 exp[f3H$ v ] 

H S v = Jl ^ (TiO-j + J2 ^ cr t a j (T kCri (25) 

«i,j>> P 

where <ii = ±1, << i,j » denotes a summation over 
nearest neighbors (diagonal bonds), and J2p denotes a 
sum over the elementary plaquettes on the square lattice. 
In the above, we have specialized to the case where spin 
couplings along horizontal and vertical bonds have been 



omitted. The critical temperature and v are given by 

e -2/3 c J 2 = ginh^J^ 

- = 2--cos- 1 (tanh(2/3 c J 2 ))- (26) 

V 7T 

We have computed Ms for this model, using system sizes 
L x L, with L = 8,12,20,40,60,80,120,200. We have 
used the Metropolis algorithm with single-spin updates, 
and up to 5.0 TO 5 sweeps over the lattice for each coupling 
constant. From this, we have extracted scaling plots of 
the peak-heights and width between peaks of M3 as a 
function of system size. From this, we directly extract 
the combinations \jv and (1 + a)/v, the final results 
for the exponents are given in Fig. [3 The agreement 




0.1 1 10 100 



a 




0.1 1 10 100 

FIG. 2: a) (1 + a) jv from FSS finite-size of M 3 for Eq. as 
a function of J2/J1. b) Same for the exponent v, computed 
directly from M3 (A) and combining results for (1 + ct)jv 
with hyperscaling (CO- c ) a as computed directly from M3 
(A) and using results for (1 + a)/u with hyperscaling (O)- 
The solid line in b) represent the analytical result Eq. 1261 
The solid lines in a) and c) are obtained from Eq. I2(il and 
using hyperscaling a = 2 — 2v. The above results thus also 
provide a check on hyperscaling in this model. 

between our simulation results and the analytical result 
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Eq. is excellent. The above demonstrates that for 
this model, the third order moment provides an excellent 
means of extracting non-universal exponents from prac- 
tically accessible system sizes, and in our view provides a 
first excellent and highly non-trivial benchmark test on 
the method. We next provide further benchmark tests 
on a number of 3D systems, before going on to extract- 
ing exponents for the 2 + 1-dimensional compact abelian 
Higgs model. 



C. Simulations of the 3D XY-, Ising-, and Ising 
(Z2) gauge-models 

We reemphasize that what follows in this section is a 
benchmark on the method of bringing out non-analytic 
thermodynamics without recourse to order parameter 
measurements, by measuring third moments of the ac- 
tion. In later sections the method will be applied to a 
lattice gauge model for which no non-local order param- 
eter currently is known, and where we will obtain precise 
values of critical exponents. In the present section, we 
will reproduce known results for the critical exponents 
a and v for the 3DXY model, and the two 3D Ising 
models, directly from third moments of the action, thus 
benchmarking the method. We will demonstrate that the 
method of using the third moment of the action is far su- 
perior to using the second moment, for practical system 
sizes accesible in numerical simulations. This is due to 
the fact that the third moment is not contaminated by 
contributions from the second moment, which unfortu- 
nately often is marred by contributions from confluent, 
or subdominant, singularities. 

The XY and Ising models are defined by the following 
two partition functions, 



Zxy — 



n d0 j 
2tt 



exp(/3i7xy), 



3=1 

H X y = ^cos(A^), 

exp(/?tfi), 

^2 ai<T i- 



Zi 
Hi 



(id) 



In addition, we consider the three-dimensional Ising lat- 
tice gauge theory [T(i| 

Zigt = exp(/3iJ IG T), 

{E/y=±l} 

HlGT = ^ UjjUjlUlkUki, 
P 

where the summation in the last expression runs over all 
plaquettes P of the lattice. We note for later use that 
the Ising model (Z 2 spin model) defined by Zi and the 



Ising gauge theory (Z 2 gauge theory) defined by Zjqt 
are dual to each other in d = 3 • 

The second and third moments M% and M3 we will 
consider for these models are given by 



Mo 



M, 



d 2 F 
d(3 2 
d 3 F 



(Gl), 
(Gl), 



(27) 



where A E (XY,I,IGT), and where I and IGT denote 
the Ising- and Ising gauge theories, respectively. In all 
these cases, scaling ~ L a l v is expected for the peak 
in the second moment, while scaling ~ T J ( 1 + a )/ L ' an d 
~ L~ x l v is expected for the third moment peak-to-peak 
height and width respectively as indicated in Fig. ^ We 
mention in passing that for the 3DXY model, we have 
(1 + a)/v = 1.467, while (1 +a)/u = 1.763 for the 3D 
Ising model and 3DZ 2 lattice gauge theory. The latter 
follows from the fact that models which are connected 
by duality transformations have identical values of a, v, 
since these two exponents can be obtained from scaling of 
the free energy. On the other hand, certain combinations 
of the remaining critical exponents, which depend on the 
degrees of freedom one chooses to describe the transition 
with, remain invariant. These invariant combinations are 
given by 0] 
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n 



2/3 + 7, PV+l). 



(28) 



The first is a consequence of Fisher's scaling law, the 
other two follow from Rushbrookes scaling law. 

Fig. [3]shows the second and third moments for the 3D 
Ising model for system sizes L — 12, 20, 32, 40 as a func- 
tion of (3, while Fig. ^shows the corresponding quantities 
for the 3DXY model. 



D. Finite-size scaling of AI2 and Mz 

We first consider the second moments of the 3DXY 
model and the 3D Ising model. These exhibit peaks at 
T c which in principle are amenable to finite size scaling. 
Fig. 03 upper panel shows the finite-size scaling plots of 
the peaks in the second moment of both models. As 
is seen from the figure, the peaks grow as L increases. 
However, it is clear that none of the scaling plots have 
reached their asymptotic behaviour. 

For the 3DXY model, the apparent scaling for small 
to intermediate values of L is clearly spurious, as it levels 
off for large system sizes. In principle, a negative slope 
should eventually be obtained for asymptotically large 
system sizes, but the general experience is that impracti- 
cally large system sizes are required to see this, let alone 
estimate the slightly negative value of a with any preci- 
sion. 

For the 3D Ising model, the situation is different in 
that the peak height grows steadily as L increases, even- 
tually approaching a straight line on a double-logarithmic 



10 




i 1 1 — 

L12 i — i- 
L20 : x 
« L32 
•x L40 r □ 
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fB 1 r 
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L20 : x : 
L32 ' x 
L40 ; :: 
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_i a i 



0.219 0.22 0.221 0.222 0.223 0.224 0.225 0.226 
P 



500 



-500 



-1000 



-1500 




0.44 0.445 0.45 0.455 0.46 0.465 0.47 

P 



FIG. 3: Second moment (upper panel) and third moment 
(lower panel) of the action for the 3D Ising model for system 
sizes L = 12, 20, 32, 40. 



FIG. 4: Second moment (upper panel) and third moment 
(lower panel) of the action for the 3DXY model for system 
sizes L = 12, 20, 32, 40. 



plot. Although the quality of the scaling improves with 
increasing L, an approximate evaluation of the slope 
yields ajv = 0.3 whereas the known value is a/v = 0.175. 
Hence, it is clear that inaccesibly large systems are re- 
quired to obtain the specific heat exponent with any ac- 
curacy from the second moment scaling analysis. 

Consider now the third moment of the action for the 
3DXY and Ising models. The results are given in Fig. 
middle panel. It is obvious that the quality of the scaling 
in both cases is vastly improved compared to the results 
obtained from the second moment. This demonstrates 
rather convincingly that finite-size scaling of the third 
moment of the action, i.e. a purely thermodynamic mea- 
surement without recourse to order parameters, suffices 
to bring out that a singular part of the free energy exists 



in both cases. This would have been difficult to conclude 
based on the second moment of the action, at least for 
the 3DXY model. Note that scaling such as this would 
not have been found had there not been a non-analytic 
part of the free energy. Thus, the quality of the scaling 
alone suffices to demonstrate the existence of a phase- 
transition. Moreover, the double-peak structure in M3 
permits a separation of a and v without recourse to hy- 
perscaling. The results for the widths between the peaks 
in the histrograms are given in Fig. [5] bottom panel. 
Note the good quality of the scaling. 

The benchmark results for the three models are given 
in Tabled We see that the error bars of these benchmark 
values all include the known values of a and v. Finally, 
the exponents obtained for the 3D Ising spin model and 
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FIG. 5: Finite-size scaling plots for the 3DXY, 3D Ising, and 
3D Zi lattice gauge model denned in Eq. 1271 a) The peaks 
in the second moment of the action demonstrating that very 
large systems sizes are necessary to obtain correct scaling re- 
sults. For clarity, two of the plots are lowered with a multi- 
plicative factor of 0.7 and 0.45 for the Ising- and Ising lattice 
gauge models, respectively, b) The peak to peak values of the 
third moment of the action, demonstrating that correct scal- 
ing results are obtained for much smaller system sizes than 
what is the case for the second moment. Note also the close 
similarity of the scaling results for the Ising and Z2 lattice 
gauge models, which should be indentical, up to constants, 
due to duality. This serves as a quality check on the simula- 
tions, c) The width between peaks in the third moment. 



the 3-D Ising gauge theory, are within error bars found 
to be the same, as they should be from duality of the 
two theories in three dimensions. Since the degrees of 
freedom used in the simulations are vastly different, this 
serves as a highly nontrivial quality check on the simula- 
tions. 



TABLE I: Critical exponents extracted for the benchmark 
models, where a is computed by combining the (1 + a)jv and 
1/v results. 



1-5 - 


Model 


(H-a)/v 


V 


a 




XY 


1.46 ±0.01 


0.67 ±0.01 


-0.01 ±0.01 


m 


Ising 


1.77 ±0.01 


0.63 ±0.01 


0.11 ±0.01 




Z2 Gau{ 


;e 1.78 ±0.02 


0.63 ±0.01 


0.12 ±0.02 



E. 



Benchmark simulations of a 3DXY model 
coupled to a Z2 clock model 



To investigate the scaling behaviour of the third mo- 
ment of the action in a well studied 3-D theory [49j with 
two coupling constants and two fixed points of type Z% 
and XY, we have considered the model 



Z 

Hxy 
Hzq 



3=1 

= ^Tcos(V^) 



exp [pH XY + hH Zq ] 



(29) 



defined on a d = 3 dimensional cubic lattice with N 
sites. The phase 6j resides on every site j and q is integer 
valued. 

Let us consider the q = 2 theory. By simple inspec- 
tion of the model we find that the limit h — > 00 leads 
to a twofold symmetry constraint on the phase, and the 
model becomes exactly the 3-D Ising model. The limits 
(3 — > 00 and (3 = are trivial theories. From renormal- 
ization group theory, the q = 2 model is found to have a 
XY phase transition with h = 0, while it exhibits Ising 
exponents for any other h ^ 0(at least in the vicinity of 
h = 0.) 

We have performed Monte Carlo simulations on the 
model Eq. (|29fl . q = 2 with system sizes up to L = 32. 
The resulting exponents from the third moment scaling 
analysis are given in Fig. |SJ We measure Z2 critical expo- 
nents down to very small values of h. However, it is clear 
that the critical coupling approaches the 3-DXY value as 
shown in the upper panel. Note the dramatic increase in 
both (l + a)/is and 1/v for h — 0.01 in the vicinity of the 
XY fixed point, where we find that the exponents become 
even less XY like. The value for a shows a more smooth 
crossover behavior from 3DXY value to 3D Ising value 
as h is increased from 0. However, the scaling appears 
not to have reached its asymptotic behaviour and the 
slopes are clearly decreasing for increasing system sizes. 
We expect this to be due to a cross-over regime in which 
the cross-over exponent contaminates our estimates, but 
that all the combinations of exponents eventually will 
converge towards Ising values, when h =/= 0. 
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FIG. 6: Upper panel: The phase diagram of the model Eq. 
1291 The open circles denote points where only the phase tran- 
sition has been located, while the filled triangles denote points 
where also the critical exponents have been measured. Lower 
panels: The three lower panels show the combination of ex- 
ponents (1 + ot)/u, as well as v and a, as a function of h//3. 
Note the rapid crossover in a from the 3DXY value at h = 
to the 3D Ising value as h is increased. The dotted and hori- 
zontal lines denote 3DXY and 3D Ising values, respectively, 
for the various quantities. 



F. Simulations of the 3D compact abelian Higgs 
model, q — 2 

We next consider the third moment of the model de- 
fined by Eq. (J5J. The model with q = 2 has previ- 
ously been studied by Monte Carlo simulations on small 
systems L 3 , with L = 2,3,4 many years ago [5(| us- 
ing second moments of the action. A short version of 



the results to be presented in this section, has already 
appeared [l9| . To our knowledge, no simulations have 
been performed in d = 3 between those presented in 
Ref. |5(J and the much more recent large-scale simula- 
tions that have been performed Due to our obser- 
vation that scaling of the second moment requires much 
larger systems than those studied in |50j | in order to ob- 
tain reliable results for exponents, a revisit to the prob- 
lem seems very appropriate. This is particularly true 
given the importance the model now has acquired in 
condensed matter system as a real model for strongly 
correlated quantum systems at zero temperature in two 
spatial dimensions 0, 0, 0, |2^| . In Ref. |2^, precisely 
such an effective model with q = 2 was derived from a 
proposed microscopic description of charge-fractionalized 
phases in strongly correlated systems. Specifically, the 
confined phase of the q = 2 model was interpreted as 
a Mott-Hubbard insulating phase, while the deconfincd- 
Higgs phase was interpreted as a charge-fractionalized 
insulating phase. 

The critical exponents obtained from third moment 
FSS analysis are presented in Fig. \7\ The exponents 
a, v, and the combination (1 + a)/v vary continuously 
along the critical line showing 3DXY and Zi universal- 
ity in the k — > oo and (3 — > oo limits, respectively. The 
Zi behaviour persists deep into the phase diagram, while 
we find a broad nonuniversal area in the large-K region. 
These two areas are joined by a peak in the exponents. 
To check the dependence of the trajectory, we performed 
simulations for several slopes at the two extremas, us- 
ing a = oo, a = 1, a = — 1 and the direction given by 
the diagonalization of the fluctuation matrix. Within 
error bars, all slopes consistently produced the same ex- 
ponents. 

The results appear to rule out that Z^- and XF-critical 
behaviors are isolated points at the extreme ends of the 
critical line. However, a feasible suggestion could be that 
two types of universality, Z2 and XY, are separated at a 
multicritical point on the critical line. We believe this to 
be ruled out by the strong deviation in (1 + a)/v from 
Zi- and XY- values at intermediate k/ (3. On balance, we 
thus conclude that the model Eq. defines a fixed-line 
theory, rather than exhibiting two scaling regimes sepa- 
rated by a multicritical point. However, the Zi character 
of the confinement-deconfinement transition persists to 
surprisingly large values of k/ (3 on the critical line, see 
also Fig. 5 of Re f. l50t Fixed-line theories in 2 + 1 dimen- 
sions are known [5l| , and non-universal exponents imply 
the existence of marginal operators in Eq. ©, yet to be 
identified. 



G. Simulations of the 3D compact abelian Higgs 
model, q — 3 

The Z q spin model and g-state Potts model are easily 
seen to be equivalent for q = 2 and q = 3 jE^. For the 
g-state Potts model, it is known that when one general- 
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FIG. 7: a) (l+a)/v from FSS finite-size of Ms for Eq. for 
q = 2. Note the variation relative to the i?2-limit 1.76 (dot- 
ted horizontal line) and the J7(l)-limit 1.467 (solid horizontal 
line), b) Same for the exponent u, computed directly from 
M3 (A) and combining results for (1 + a)/v with hyperscal- 
ing (O)- c ) Q as computed directly from M3 (A) and using 
results for (1 + a)jv with hyperscaling (O)- The maximum 
and minimum in a) have been obtained by crossing the critical 
line along the trajectory (3(n) = /3 C + a (k — n c ) with a — 00 
(A), o = l (■), and a = —1 (A) using C = 0.665, k c = 2.125 
(max.), and /3 C = 0.525, k c — 5.0 (min.). 



izes q to be real- valued, the phase transition in the model 
in d = 3 changes from continuous to discontinuous when 
q is increased beyond the value q = 2.625 |53|. Hence, 
both the lattice gauge theory and the Z3 spin model 
have first order phase transitions. In the limits (3 — > 00 
and k — > 00 the model defined by @ reduces to the Z3- 
spin model and 3DXY model respectively. If these phase 
transitions survive on the critical line for finite coupling 
values, a tricritical point joining the first and second or- 
der critical lines is expected to exist. As we shall see, this 
is precisely what happens. 

Therefore, we reach the important conclusion that the 
model exhibits a first order phase-transition not only for 
(3 = 00, Eq. (0), but also for finite values of (3. To iden- 
tify and investigate the first order phase transitions we 
have computed histograms of the action at critical 
coupling constants. These have a double-peak structure 
due to the two coexistiting phases that characterize a 
first-order phase transition. The histograms have been 
produced from large scale Monte-Carlo simulations fol- 
lowed by muliti-histogram reweighting and jackknife er- 
ror analysis The results are shown in Fig. |HJ 




<S>/N 

FIG. 8: Normalized histograms of the action @ for q = 3 
at the first order phase transition with various critical cou- 
plings and system sizes L = 8, 12, 16, 20, 24, 32, 40, 64. The 
histogram height increases with system size. The horizontal 
axis shows the average value of the action per unit volume, 
and the histograms have been shifted horizontally for clarity. 

A first-order phase transition is characterized by two 
coexisting phases with the same free energy. There 
should therefore exist domain walls separating the two 
phases. The area of the domain walls is related to the 
energy difference AF(L) required to keep the two phases 
separated through the expression 

AF(L) = ]nP(S, L) max - lnP(5, L) min ~ L d -\ (30) 

where P(S, L) is the probability for a value S of the ac- 
tion in a system of size L d , and L d ~ x is the cross-section 
area between the ordered and disordered phases j5^|. The 
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results shown in FigO confirm this for system sizes L > 
24. 

The discontinuity of the action or equivalently the 
width between the peaks in the histograms is, strictly 
speaking, not the latent heat. Nonetheless, it can be 
taken as a reasonably accurate measure of this quantity 
along a sufficiently small part of the phase transtion line. 
It is important to ascertain whether or not the hrst-order 
character of systems persists in the thermodynamic limit. 
Hence, we plot the discontinuity of the action as a func- 
tion of system size in Fig. [5J upper panel, and find that 
the discontinuity of the action is finite at least up to 
k c = 1.65. The first order phase transition becomes in- 
creasingly weaker, i.e. the latent heat in the transition is 
reduced as we approach k c — 1.65 from below. 



dynamic limit of this quantity, and plotted them as a 
function of the critical coupling in Fig. ^| A linear 
extrapolation yields an estimate for the tricritical point 
fttri/Ari = 1.39 ± 0.06 corresponding to (0trU ^tri) = 
(1.23 ±0,03, 1.73 ±0.03). 




50 60 70 



FIG. 9: Upper panel: The discontinuity of the action J5J as 
a function of system size. Note that the plots seem to be 
converging towards a finite energy value. Lower panel: The 
energy difference between the two coexisting phases AF(L) 
as a function of system size, same labeling as in the upper 
panel. 

Approaching the tricritical point along the first or- 
der line, the discontinuity in the action, equivalently the 
width between the peaks in the histograms, must van- 
ish. From Fig. [5] we have extrapolated the thermo- 




FIG. 10: Width of peaks in histograms as a function of k c / (3 C . 
The solid line is a linear fit to the data points for k c //3 c > 0.8. 

On the other side of the tricritical point we have per- 
formed a third moment analysis. The exponents are ex- 
pected to be those of 3DXY in the large- k limit. Strictly 
speaking, critical exponents are properties of second or- 
der phase transitions. For a first order transition the 
correlation length stays finite and hence there is no 
scale-invariance. However, bounds for some of the ex- 
ponents can be obtained and the limits of these bounds 
correspond to the exponents one would get by relax- 
ing the definitions and formally considering first order 
transitions ,43 . One finds that the limits are a = 1 and 
v = 1/d respectivly. Hence when approaching the tricrit- 
ical point we expect (1 + a)/i> = 6 and v = 1/3. 

The critical exponents, given in Fig. 1111 agree with 
the expectations. The exponents are XY like in the large 
k limit while in the region 2.4 < k < 50 the quantity 
(1 + a)/v is lower than the XY value and hence exhibits 
nonuniversiality similar to what we found for q — 2. At 
approximately n ~ 2.4 the exponents rise abruptly to- 
wards the first order limiting values. Due to super-critical 
slowing down, obtaining good quality third moment scal- 
ing plots becomes difficult when approaching the tricrit- 
ical point. Thus, locating the critical point by deciding 
where the exponents have reached the values one expect 
in the limit of a first order phase transition, is not feasi- 
ble. 

For q = 3, which is a special case in this context (see 
below), this is very much like what is now known to hap- 
pen in the non-compact abclian Higgs model in d = 3. 
In that system, a first order phase transition character- 
istic of type-I superconductivity at small values of the 
Ginzburg-Landau parameter, is converted to a second 
order phase transition characteristic of type-II supercon- 
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FIG. 11: a) (1 + a)/v from FSS finite-size of M 3 for Eq. 
@ for q = 3. Note the variation relative to the !7(l)-limit 
1.467 (solid horizontal line) and the violent behaviour near 
the critical point n tr i- b) Same for the exponent v, computed 
directly from M3 (A) and combining results for(l + a)/^ with 
hyperscaling (CO- c ) a as computed directly from M3 (A) and 
combining results for (1 + a)jv with hyperscaling (CO- 



ductivity. The tricritical value of the Ginzburg-Landau 
parameter where this change in the character of the phase 
transition occurs, has recently been determined with pre- 
cision in large-scale Monte Carlo simulations, to be given 
by 0.8/v^ m- This 

is in remarkable agreement with 
previous analytical results using duality arguments |56j| . 



H. The phase diagram for q — 2, 3, 4 

Let us summarize what has been discussed above. The 
phase structure of the model is given in Fig. From 
the early Monte-Carlo simulations on the model [5(j for 
q = 2, it is known that the critical line (3 c (k) approaches 
linif^oo [3 C = 0.454 while when (3 — > 00, there is a critical 
value of k given by k c = 0.761. As q increases, the verti- 
cal part of the phase-transition line moves up in k, while 
for the values of q we have considered, the lines are criti- 
cal except for the case q = 3 which is a first order line for 
large values of /?. For any q, however, we know that the 
limit k — ► 00 represents the U(l) limit, which must ex- 
hibit a second order phase transition. Hence, the line for 
q = 3 in Fig. 1 121 must contain a tricritical point. The first 
order line for large (3 terminates at a tricritical point, and 
is second order for larger values of n. For q — 4, which 
is also shown, the phase-transition line is second order. 
The limit (3 — > 00 corresponds to the d = 3 Z4 lattice 
model which is dual to the Z4 spin model. Using the fact 
that the symmetry group Z4 = Z2 <8 Z2 , the universality 
class of the q = 4 phase transition line is expected to 
interpolate between the Ising case in the limit (3 — ► 00 
and the 3dXY universality class in the limit k — > 00. 




FIG. 12: The phase diagram for the d = 3 compact abelian 
Higgs model in three dimensions for q = 2,3, 4, 5. All lines are 
critical for all values of k except for the case q = 3, which is 
first order for k < K t ri and second order otherwise. The thick 
solid portion of the q = 3-line denotes a first-order transition. 
Detailed finite-size scaling analysis of the line q = 2 shows that 
the critical exponents a and v vary continuously along the 
critical line. This makes the theory defined by Eq. a fixed- 
line theory with non-universal critical exponents, as opposed 
to a fixed-point theory with universal exponents. The critical 
line approaches the 3DXY value /3 C = 0.453 as k — > 00 for all 
values of q £ N. 

Based on what is presented above, we conclude that 
for q 7^ 3, the model Eq. J5J exhibits a critical line (3 c (n) 
for all values of k, while for q = 3 the phase transition 
of the model is first order for small values of k and is 
converted to a second order phase transition when k is 
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increased beyond a finite tricritical value. 

V. CONCLUSIONS 

In this paper we have, in addition to a number of 
bench-mark models, considered various lattice gauge 
models in 2 + 1 Euclidean dimensions, with particular 
emphasis on the compact abelian Higgs model with inte- 
ger gauge charge q > 2. We have seen that that for all 
1 > 2; q ^ 3, the phase transition for the model in Eq. 
P)l is second order such that the entire phase transition 
line /3pt(k) is critical, i.e. /3pt(^) = f3 c (n). The expo- 
nents a and v which are given by the 3DXY values in 
the limit k» f3, and global Z q values in the limit [3 n. 
For the special case q = 2 and for intermediate values of 
(3 and k on the critical line, the exponents vary contin- 
uously from the value (1 + a)/v = 1.76 in the Z 2 -limit 
j3 > k to the value (1 +a)/v = 1.467 in the U(l) limit 
s » |3. This constitutes a rare example of a fixed-line 
theory with non-universal exponents in a three dimen- 
sional system [5j|. What the connection to the work of 
Refs. 113 is, where a Kosterlitz-Thouless like phase 
transition of unbinding of monopoles in a three dimen- 
sional compact U(l) gauge theory was found, remains to 
be investigated. More work is also required to elucidate 
the special role of q = 3 in the compact abelian Higgs 
lattice gauge model. Note that this is fundamentally dif- 
ferent from what happens in the g-state Potts model, to 
which the Z q spin model is equivalent only for q — 2,3. 
In constrast to the compact abelian Higgs model, the q- 
state Potts model has a first order phase transition for 
all q > 2.625 [j^. 

Based on the connection between the model Eq. @ 
for q = 2 and a recently proposed model of fractional- 
ized phases in strongly correlated systems j2^| which is 
essentially given by Eq. |J2J, we propose that the uni- 
versality class of the putative quantum phase transition 
from a Mott-Hubbard insulator to a charge- fractionalized 
insulator which the model is supposed to describe, is in 
the universality class of the q = 2 compact abelian Higgs 
model characterized by a fixed line of nonuniversal crit- 



ical exponents varying continuously between the values 
(1 + a)/u = 1.76 in the limit (3 — > oo, k finite, and 
(1 + a)/v = 1.467 in the limit k —> oo, (3 finite. How- 
ever, we have found that over a significant portion of the 
critical line of the q — 2 compact abelian Higgs model, 
the exponents take on values consistent with those of the 
3D Ising model. If microscopic models describing such 
Mott insulator- fractionalized insulator transitions yield 
the q — 2 compact abelian Higgs model with sufficiently 
large values of (3 and small values of n as an effective the- 
ory, then the resulting insulator-insulator transition is in 
the 3D Ising universality class |28| . 

It appears that the study of matter-coupled compact 
gauge theories would benefit greatly from a precise for- 
mulation of a non-local order parameter represented by 
generalized rigidity for such models, as a substitute for 
the Wilson- or Polyakov loops which have proved use- 
ful in the pure gauge theories, but which are always 
perimeter-law bounded in the matter coupled case when 
the symmetry group of the matter field is a subgroup of 
the symmetry group for the gauge field. 
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